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Abstract 

We present a detailed study of the effects of the initial distribution on the kinetic evolution of the 
irreversible reaction A + B — > in one dimension. Our analytic as well as numerical work is based on 
a reaction-diffusion model of this reaction. We focus on the role of initial density fluctuations in the 
creation of the macroscopic patterns that lead to the well-known kinetic anomalies in this system. In 
particular, we discuss the role of the long wavelength components of the initial fluctuations in determining 
the long-time behavior of the system. We note that the frequently studied random initial distribution is 
but one of a variety of possible distributions leading to interesting anomalous behavior. Our discussion 
includes an initial distribution with correlated A — B pairs and one in which the initial distribution forms 
a fractal pattern. The former is an example of a distribution whose long wavelength components are 
suppressed, while the latter exemplifies one whose long wavelength components are enhanced, relative to 
those of the random distribution, 

1 Introduction 

Diffusion-limited reactions in low-dimensional geometries are well known to exhibit "anomalous" kinetics 
different from those predicted by the standard law of mass action [Ovchinnikov & Zeldovich, 1978; Toussaint 
& Wilzcek, 1983; Kang & Redner, 1984; Kopelman, 1988; Burlatskii et at, 1989]. An example of anomalous 
kinetics is exhibited by the diffusion-limited irreversible reaction A + B — ► 0. If the distribution of each 
reactant is spatially homogeneous, the rate laws for the global reactant densities pA(t) and PB(t) are pa = 
pB = —KpaPb, where if is a constant global rate coefficient. If pa(0) — Pb{0) = p{0) then the densities 
of the two species are equal at all times, pji(t) = PB{t) = p{t), and we can write p = —Kp 2 . In integrated 
form p~ l {t) = /o(0) _1 + Kt; at long times p ~ (iff) -1 . We call this behavior "classical." However, the 
actual asymptotic rate law in an infinite volume in dimensions d < 4 for an initially random distribution of 
reactants is instead p = — Kp^ 1+i / d ^ [Ovchinnikov & Zeldovich, 1978; Toussaint & Wilzcek, 1983; Kang & 
Redner, 1984; Kopelman, 1988; Burlatskii et at, 1989]. The associated asymptotic time dependence of the 
density is then p ~ £~ d / 4 (with logarithmic corrections for d = 4). 

The physical origin of the "anomalies" is in general clearly understood. The usual law of mass action 
assumes a spatially homogeneous mixture of reactants; however, in low dimensions diffusion is not an efficient 
mixing mechanism and may be incapable of smoothing out reactant concentration fluctuations that are 
created, for example, by nonhomogeneous initial conditions, by external reactant sources if the system is 
open, or by the reaction itself. 

Consider the A + B — > reaction starting from an initially random distribution of A and B. Such 
a distribution exhibits local density fluctuations, that is, certain regions of the system are relatively rich 
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in species A while others are relatively rich in species B. By "efficient" mixing we mean mixing that is 
rapid relative to the local reaction rate. Thus, if efficient mixing (by whatever mechanism) occurs, such 
local fluctuations are rapidly suppressed. If mixing is not efficient as, for example, in a diffusion-limited 
reaction in low dimensions, then the reaction will cause the local minority species to be eliminated extremely 
rapidly by the local majority species, and diffusion is not able to homogenize the system. There follows an 
evolution of regions in which essentially only one or the other of the species is present, and any attempt of 
the minority species to diffuse into this region leads to its rapid disappearance by reaction. The reaction then 
only takes place at the interfaces of aggregates of different species. As time proceeds (and the concentrations 
of reactants decrease), the regions in space that are occupied by reactants of essentially only one species 
grow in size, and the total interface area at which the reaction can occur decreases. This leads to ever more 
effective spatial segregation of the species and an associated slowing down of the reaction relative to its rate 
in a homogeneous mixture. The system goes through a hierarchy of anomalous kinetic behaviors [Argyrakis 
et al, f993; Lindenberg et al, 1996] to arrive asymptotically at the so-called Zeldovich regime [Ovchinnikov 
& Zeldovich, 1978; Burlatskii et al, 1989] characterized by the time dependence p ~ t~ d / 4 for d < 4. 

This discussion punctuates three interesting features of these irreversible diffusion-limited reactions in 
constrained (low dimensional) geometries. One is that the initial distribution determines the evolution of the 
system for all time: the future of the reaction is imprinted in the spatial distribution at one given ( "initial" ) 
moment. A second is the fact that it is necessary to start with disorder in order to produce an ordered state 
at a later time: the formation of ever larger aggregates of a single species is a direct consequence of initial 
fluctuations that are not homogenized. A third is the realization that a random distribution is but one of 
many possible initial distributions, and it is certainly not the most easily realized in practice. The historical 
focus on the random distribution has perhaps obscured the fact that other initial distributions lead to other 
behaviors, some "less anomalous" and others even "more anomalous" than the Zeldovich law p ~ i _d / 4 in 
their deviation from classical kinetics. It is this last point that we focus on in this paper. 

It is reasonable to expect that the suppression or enhancement of initial density fluctuations will affect 
the kinetic progression and change even the asymptotic behavior of the A + B — > system. In particular, an 
initial distribution that is more homogeneous or mixed than a random one is expected to lead to behavior 
closer to classical; conversely, an initial distribution in which fluctuations are more prominent is expected 
to lead to greater deviations from the law of mass action. Our work expands on previous concern with 
this problem by a number of authors [Vitukhnovsky et al, 1988; Abramson et al., 1994; Abramson, 1995; 
Lindenberg et al, 1994, 1996; Sancho et al, 1996]. 

A "mixing mechanism" that we must mention at the outset because it necessarily affects the numerical 
simulations that parallel our theoretical analysis (and because it occurs in physical systems) is provided 
by the boundaries of a finite system of volume V — L d , where L is the linear dimension of the system. 
Regardless of the precise nature of the boundaries (reflective, periodic), segregation of species is in any case 
limited by the finite size of the system. With periodic boundaries a reactant "perceives" this finite size at a 
time [Lindenberg et al, 1996] 

t f ~L 2 /8n 2 D, (1) 

where D is the diffusion coefficient for that reactant (the numerical factors are slightly different for other 
boundary conditions ).[] 

In Sec. H we summarize the model and our general theoretical framework. In Sec. [| we study various 
initial conditions and present for each case our theoretical and numerical results. A summary of conclusions 
is presented in Sec. |[ 

2 Theoretical Framework 

We base our analysis on the reaction-diffusion equations utilized in many theoretical studies of this problem. 
It is important to point out that our numerical results are also based on these same equations. This differs 
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from the more usual, more microscopic, approach of simulating the system as a lattice gas in which the 
reactants are random walkers and in which the reaction is an annihilation that occurs with certainty or with 
some probability when reactants step on the same lattice site. Simulations at this level of detail require large 
computer capabilities and long times to reach asymptotic behaviors [Argyrakis et al, 1993; Lindenberg et 
al., 1994, 1996]. Reaction-diffusion equations represent mesoscopic approximations to this more microscopic 
picture; such approximations can therefore only provide a theoretical backdrop over spatial and temporal 
ranges that are not too detailed and for densities that are not too low. The results of direct numerical solution 
of the reaction-diffusion equations are of course constrained in the same way as are the equations themselves. 
On the other hand, such solutions, being more mesoscopic from the outset, require less computational 
intensity and length to cover the hierarchies of behaviors and to achieve asymptotic behaviors. Although 
this approach of course provides no information on the limits of validity of reaction-diffusion models, they 
do provide a simpler and less costly backdrop against which to check the theoretical analysis of the model. 
The reaction-diffusion equations for the local densities pA(x,t) and pg(x,t) of the two species are: 

= DV p A - KpAPB 

at 

^ = DV 2 p B -Kp A p B (2) 
The reaction takes place under stoichiometric conditions for the global densities, 

PA(t) = J dx p A {x, t) = p B {t) = p(t), (3) 

so that the densities of both species remain equal forever. D is the diffusion coefficient, which we assume 
to be the same for both species, K is the local reaction rate constant, and V is the system volume. Note 
that we are using the same symbol, K, for the local rate coefficient as we did for the global rate coefficient. 
When the global rate law is classical, the two are in fact the same; when it is not, the global rate is in 
any case not quadratic in the global density and the "rate coefficient" in a quadratic rate law is then not 
even a constant. We vary the values of the parameters in our analysis, but recognize that (in appropriate 
dimensionless units) a diffusion-limited reaction occurs if D/K -C 1 - indeed, a strictly diffusion-limited 
reaction requires K —> oo. With K finite (no matter how large), when the density is sufficiently low the 
reaction becomes reaction- rather than diffusion-limited 

A few comments concerning notation are useful at this point. The local densities pi(x,t) evolve from 
initial configurations Pi(x, 0). In some cases these initial configurations are chosen from a distribution of 
different configurations (e.g., a random distribution) and it might be appropriate to perform an average 
over the initial distribution [Lindenberg et al., 1996]. In other cases the initial configuration may be a very 
specific one, e.g. as in the study of the evolution of fronts [Koo & Kopelman, 1991]. We denote an average 
of a quantity f(x,t) over the initial distribution by a bracket: (f(x,t)). If we are dealing with a specific 
initial configuration then the initial distribution is a delta function and a quantity "averaged over the initial 
distribution" is simply the same as the unaveraged quantity, (f(x,t)) = f(x,t). 

Global quantities such as the global density are related to local ones by the volume average indicated in 
Eq. §). We also use a single bracket notation to denote an average over volume. The use of a single bracket 
to denote one or the other or both averages should not create confusion since the context makes it clear 
what averages have been taken. In any case, a dependence only on the time indicates that a volume average 
has been performed, fit) = (f(x,t)). 

The first impetus toward finding the global density from the reaction-diffusion equation might be to 
average Eq. (|^) over initial distribution and volume so as to obtain an equation for the rate of change of the 
global density. However, this of course introduces a two-density average on the right hand side: 

d(p A {x,t)) d(p B (x,t)) 

dt = dt = ~ K {PAPB) ' (4) 

where we have noted that V 2 p(i) = since p(t) is independent of x. One might then write an equation 
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for the two-density product and again perform an average, which introduces a three-density term. In many 
problems it is possible to truncate such a hierarchy and thereby close the set of equations. Such closure 
is not possible here because it would miss the many-body nature of the problem: indeed, the evolution 
of macroscopic structures in these systems would seem to indicate that no "finite body approximation" is 
appropriate, and that one needs to capture at least some aspects of the many-body effects to all orders to 
obtain a reasonable approximation. 

We note that classical behavior is captured by Eq. (Q) only if the average of the product can be broken 
into the product of the averages, (paPb) = (pa) (pb)- This separation is in general not correct and, in 
particular, not when there is species segregation. Indeed, if the species are strictly segregated then the right 
hand side is zero because it is not possible to find A and B at the same location, and thus the rate of change 
of the global densities vanishes, as it should. If there is some species overlap, then the decay rate does not 
vanish but is slow (slower than classical). If the species are randomly or homogeneously distributed, then 
the average of the product does reduce to the product of the averages and the kinetics are classical, with a 
consequent inverse time decay of the global densities. 

A number of procedures have been introduced to deal with the many-body aspects of the problem 
[Vitukhnovsky et ai, 1988; Burlatskii et ai, 1989]. In our analytic studies we follow [Vitukhnovsky et 
ai, 1988; [Lindenberg et al., 1988, 1990; Sokolov & Blumen, 1991; Sancho et ai, 1996] by introducing the 
difference and sum variables, 

7(5 , t) = p A (x,t) - PB (x,t) ^ p% t) = p A (x,t) + PB (x,t)^ (5) 
so that the pair of equations (^) transforms into 

P) 

The difference variable "f(x,t) captures the variations and fluctuations in the spatial distribution of the 
species. If the system is thoroughly mixed, this variable is everywhere small; if the species are segregated, 
then the variable is positive in regions where A predominates and negative where B predominates; its 
variability contains information on the sizes of such regions. If the species were oppositely charged, for 
example, then this variable would represent the local net charge. Fortunately, Eq. (^) is a simple linear 
diffusion equation and can be solved exactly [Vitukhnovsky et ai, 1988 Lindenberg et ai, 1996] (and it is 
through this exact solution that one captures important many-body effects to all orders): 

j(x,t) = dx'e- Dtk2 e^- s 'h(£',0) (8) 

k 

where k = \k\. For periodic boundary conditions k = 2~Kn/L and ft is a d-tuple of integers. In the limit 
V — > oo one can transform the sum to an integral: 

l{x,t) ^^L^JdkJ dx'e- D ^e i ^- s 'h(S',0). (9) 

The integration over k can easily be done explicitly in any dimension; in this paper we deal with one- 
dimensional systems, for which one readily obtains 

7 (x, t) = -J= J dx'e-^-^/^ix', 0), (10) 
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where we have dropped the vectorial notation for x and x' but note that they can be positive or negative. 
We return to the difference variable momentarily. 

Consider next Eq. (^) for the sum variable. This variable is the density without consideration of species 
identity, and it is in general expected to be a smoother function than the difference variable. Let us average 
Eq. (0) over the volume and over the initial distribution of reactants: 

d M = - K m ( 7 2)] . (ii) 

The entire subsequent analysis of the kinetic behavior of the average global density p(t) = (p(x,t)) is based 
on this equation, and, in particular, on the behavior of the function 

G(t)^( 7 2 (x,t)) (12) 

(which one knows exactly), and on the relation between (p 2 (x, t)) and p(t). We shall call G{t) the segre- 
gation function, since the difference variable (and consequently this function) vanishes for a homogeneous 
distribution. 

As we demonstrate explicitly below, the decay of G(t) with time (through diffusion) is not simply of a 
fixed inverse power form over the entire history of the process. Furthermore, in general there is no simple 
relation that holds for the entire progression of the reaction between the average of the square of the local 
density and the average density itself. Consequently, there is no exact kinetic description of the system 
in terms of only the global density: indeed, Eq. ( [TT|) is again but the first member of an infinite moment 
hierarchy. Nevertheless, one attempts to extract a kinetic description from Eq. (^), at least approximately, 
recognizing that in any case a single kinetic description does not hold for all times. It has therefore become 
common practice to speak of kinetic behaviors that are valid over some substantial time range, and crossovers 
from one type of behavior to another. If the "crossover" regimes are short compared to the regimes where one 
type of behavior persists, one can then characterize the system in terms of a sequence of kinetic behaviors and 
transitions between these behaviors. Such a description turns out to be possible in many cases. Toward this 
purpose, in [Lindenberg et al., 1996] we presented an extensive discussion concerning the relation between 
the average of the square of the density, (p 2 (x, i)), and the average density, p(t) = (p(x,t)). In particular, 
we noted that in general (p 2 (:?,£)) ^ p 2 (t)- For a random distribution of reactants or of reactant pairs 

(p 2 {x,t)) = p 2 (t) + p max p(t), (13) 

where p m ax is the "maximum density" [that is, p(t) < p max ] specified later for our particular reaction- 
diffusion approach^ On the other hand, for a segregated distribution as occurs in the Zeldovich regime, and 
also for a homogeneous distribution (these two "opposite" cases share this second moment property) 

(p 2 (x,t))=p 2 (t). (14) 

In any case, the balance of the three terms in Eq. ( |ll| ) determines the behavior of the global density. A 
general point that we return to subsequently in the context of particular cases can already be made here: 
the presence of the "driver" G(t) = (j 2 ) causes deviations from classical behavior; as already mentioned, 
a homogeneous distribution of reactants at time t is associated with a vanishing G(t) and with classical 
kinetics. 

The behavior of G(t) for all times and, in particular, how rapidly it vanishes at long times, is completely 
determined by the initial configuration of reactants. It is convenient to write 

GW^E^V 2 ^' 2 (15) 

k 



2 Prom a molecular viewpoint pmax is the inverse of the volume of one molecule. In a lattice approach it is the maximum 
number of walkers that can simultaneously occupy each lattice site. 
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or, in the large volume limit, 



G(t) = j^yl J dkS(k)e~ 2Dtk \ (16) 

where we have used Eq. @. In either case S(k) is the Fourier transform of the initial correlation function, 

S(k) = J dxe~ fk - s ( 7 (x, 0) 7 (0, 0)) . (17) 

S(k) is called the structure function of the initial configuration. For t > t a , where 

t a =(Dk 2 a )-\ (18) 

G(t) is essentially determined by the behavior of the structure function for k < k a . In other words, the long- 
time behavior of the segregation function defined by the difference variable is dictated by the long-wavelength 
components of the initial distribution (and also note that this behavior is entirely determined by diffusion 
and is not influenced by the reaction) . In this paper we consider a variety of initial conditions characterized 
by different forms of the structure function. Familiar cases include the random initial distribution in an 
infinite volume, for which S(k) = const for k > 0, and a random initial distribution of A~ B pairs separated 
by a distance c in an infinite volume, for which S(k) ~ k 2 for k < c" 1 . We return to these cases subsequently. 
In general, if S(k) ~ k a for k < k a , then a rescaling of variables in Eq. ( |l6| ) immediately shows that for 
t > t a 

G{t)~t-^. (19) 

Note that a may be positive or negative (or zero). If a is sufficiently large (indicating a drastic cutoff of 
the long wavelength components), specifically if a > 4 — d, then G(t) decays more rapidly than t~ 2 and it 
rapidly ceases to play a role in the kinetics. The behavior then becomes classical [Vitukhnovsky et ai, 1988; 
Sancho et at, 1996]. Conversely, any initial distribution that emphasizes the long wavelength components 
in the density difference leads to a slower decay of the "driver" (a < 0) and hence to greater deviations 
from classical behavior. We present examples of both cases. Note that it would in any case at this point be 
premature to conclude from Eq. ( |l9| ) that p(t) decays as £-( Q + d )/ 4 - at the very least, the moment relation 
( |l4|) must be established first. 

Ultimately we are constrained by the finite systems in which our numerical work is carried out; this 
means that eventually G(t) vanishes exponentially, reflecting the fact that there are no contributions to it 
whose wavelengths are longer than the size of the system. 

The results of our analysis of the reaction-diffusion equations will be compared with the numerical 
integration of Eq. (|J). This integration is carried out by discretizing the volume V = L into N cells of 
size Ax (L — NAx) and time into intervals At. The Laplace operator is represented in the standard 
symmetrical form with respect to the mesh size Ax, and the time integration follows a standard first order 
Euler algorithm. Note that adjustment of the discretization parameters Ax and At can alternatively be 
accomplished by adjusting the parameters D, K, and N. We always set K ^> D in order to be in the 
diffusion-limited regime for a meaningful interval of time, but otherwise set these parameters so as to be 
able to clearly illustrate the behaviors under consideration. Throughout we set Ax = 1, and At = 0.01 
except for early time calculations where on occasion we use smaller values of At for the first few hundred 
integration steps. Initially an amount po of reactant A is placed in some cells and not in others, according to 
a specified distribution. Similarly for reactant B. If each reactant is initially placed in n < N cells, then the 
initial global density is Pinit = npg/N. Note that the maximum initial density possible is then Pmax — Po> 
which occurs if every cell is initially occupied by both reactants. An average over different initial conditions 
is taken if necessary. 
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3 A Variety of Initial Conditions 



The initial pattern of fluctuations determines for all time the form of the driver G(t) in Eq. (jll]) with (|l2|). 
We establish below through specific examples that this driver decays more slowly (rapidly) when the long 
wavelength components of the initial difference fluctuations are enhanced (suppressed) , and that the kinetics 
of the A + B — > reaction are more anomalous, i.e., further from classical behavior, when the driver decays 
more slowly. When the decay of G(t) is more rapid than t~ 2 then its contribution becomes irrelevant and 
the rate law approaches the classical form. 

In this section we present the kinetic progressions for the A + B — > reaction for a number of initial 
conditions. The case of the random initial condition has been extensively studied in the literature and 
is reviewed here in detail because it serves as a benchmark for the discussion of other initial conditions 
[Argyrakis et al, 1993; Lindenberg et al, 1996]. This review allows us to establish some limitations imposed 
by the reaction-diffusion approach [Lindenberg et al, 1996]. We then present an example of an initial 
condition in which long wavelength components of the density difference fluctuations are enhanced relative 
to those of a random distribution. This initial condition, in the form of a fractal pattern, leads to more 
pronounced deviations from classical behavior. We also present an example of an initial distribution whose 
long wavelength components are more constrained than the random, namely, the case of correlated pairs 
[Lindenberg et al, 1994, 1996, Sancho et al, 1996]. 



3.1 Random initial distribution 

The random initial distribution of reactants has been the most thoroughly studied, even though it is difficult 
to implement in practice. Nevertheless, it is important because it is the distribution that was first used 
to quantify the fact that diffusion in low dimensions is not an effective mixing mechanism. It was the 
distribution that clearly showed that initial (microscopic) fluctuations are a necessary element to achieve 
(macroscopic) segregated patterns, that is, that initial inhomogeneities, with no further explicit separating 
mechanism or interactions other than the reaction itself, lead to macroscopic species segregation in low 
dimensions [Ovchinnikov & Zeldovich, 1978; Toussaint & Wilzcek, 1983; Kang & Redner, 1984; Kopclman, 
1988; Burlatskii et al, 1989]. 

We have argued elsewhere [Argyrakis et al, 1993; Lindenberg et al, 1996] that at extremely early times 
the kinetics of the system is classical, that is, that with a random distribution of reactants the global density 
obeys the classical rate law dp(t)/dt = —Kp 2 (t). We also noted that this behavior in low dimensions persists 
only for a time of order K~ x ; on this time scale the distribution already begins to deviate from the random 
distribution because of the reaction of physically contiguous A and B molecules. Because here we compare 
our analytic results with those obtained numerically according to the discretization scheme described earlier, 
we are constrained from consideration of initial concentration fluctuations on a spatial scale smaller than Ax. 
In our discretized reaction-diffusion approach A and B densities that coexist within one cell are consumed 
on the time scale K~ l , and each cell is, after a time K~ x , either empty or occupied essentially only by A or 
only by B [Sancho et al, 1996]. 

Beyond this very early time behavior, the results that we detail below are captured in Fig. [l]. We first 
describe the expected behavior and then compare it to that seen in the figure. As noted earlier, an initially 
random distribution in an infinite volume is characterized by the structure function S(k) — const for k > 0. 
In our discretized integrations over a finite volume, a random distribution is characterized by a constant 
structure function only in the range k^ < k < ko, where k^ 1 = Ax is the limit imposed by the discretization 
and Ul is the limit imposed by the finite system size. It follows from Eq. ( |Tg| ) that for times t > to, where 

(») 



the segregation function G{t) decays as t d l 2 [cf. Eq. (119)]; in fact, exact integration of Eq. ( fig) leads to 
(we return to finite volume effects in a moment) 

G(t) = Q{Dt)- d ' 2 , (21) 

where Q is a constant proportional to the initial density. 
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Figure 1: Decay of the global density in one dimension from an initially random distribution of reactants. 
In both figures TV = 512 and K = 10. (a) D = 2.0, (b) D = 0.5. The curves are the numerical solutions of 
the reaction-diffusion equations. Straight lines in (a): £ -1 / 4 and t~ x . Straight line in (b): i -1 / 4 . 



The consequent behavior of the global density predicted by Eq. (|11|) then depends on the relation between 
the average of the square of the density and the global density itself. We have noted that for a random 
distribution Eq. ( |l3| ) holds. This relation is nearly satisfied for a nearly random distribution; that is, at 
early times before the reactant distribution deviates too much from randomness the contribution linear in 
the density dominates and the rate law is of the form 



dp{t) 
dt 



-K[p 

max 

p(t) - Q(Dt) 



-d/2] 



(22) 



For d = 1 the dominant terms are those on the right hand side, and the global density thus decays as 

p(t) ~t-i/2. (23) 

This decay is characteristic of the presence of a "depletion zone" around each reactant (and, for the same 
reasons, is also characteristic for the A + A — > reaction). The critical dimension where all three terms are 
of the same order is d = 2. Note that consistency in one dimension requires that the difference between the 
two contributions on the right decay in time as t~ 3 / 2 so as to be balanced by pit). 

Deviations from randomness and, in particular, spatial segregation of reactants, begin to set in very 
quickly for d < 4 (for a detailed discussion of crossover times see [Argyrakis et at, 1993; Lindenberg et ai, 
1996]). The moment relation changes to (|l4|), and the rate law becomes 



dp(t) 
dt 



-K[p 2 (t)-Qr d l 2 ] 



(24) 



For d < 4 the dominant terms are again those on the right hand side, and the global density in one dimension 
now decays as 



P((t) 



t 



-1/4 



(25) 



("Zeldovich" behavior). The critical dimension where all three terms are of the same order is d = 4. In 



one dimension the difference between the dominant terms, and aiso pit), now decay as t -5 / 4 . The Zeldovich 
behavior reflects the well-known segregation of the species - even as the densities decrease, the sizes of the 
regions essentially occupied by only one species increase. Note that the fact that ( |l4| ) rather than (|l^) holds 
in this segregated regime indicates that the spatial distribution of the reactants within the single-species 
aggregates is not random. If it were random, a relation such as ( |l3| ) (with slightly modified coefficients) 
would hold instead. 

Eventually, when t — > tf [cf. Eq. (jjj)], the above analysis must be adjusted to reflect finite system size 
effects. The integral form Eq. ( fl6| ) can no longer be used, and one must return to the sum form Eq. (|l5|). 
The segregation function for t 3> tf decays exponentially. There is a rapid drop in the global density due 
to the increased reaction rate caused by the mixing of species forced by the boundaries of the system. A 
"homogeneous" regime is then reached where the last term in Eq. (|i"T| ) is negligible, and the rate law returns 
to the classical form 

tf-~-K P >(t), (26) 

with the associated decay 

p(t)~iKt)-\ (27) 

Figure [j] shows the evolution of the global density for a one-dimensional system with a random initial 
distribution of cells containing A and a random initial distribution of cells containing B. We generate each of 
these distributions separately but with the same number of occupied cells for each species, and then simply 
superpose them. This produces some "premixed" cells whose density decays very rapidly [Sancho et al., 
1996]. Extensive premixing is avoided by starting with a relatively low initial density. We calculate the 
resulting evolution for two different values of the diffusion coefficient. Figure |l|(a) is for a larger diffusion 
coefficient than Fig. |l](b). In both, the Zeldovich decay regime is clearly seen over a considerable portion of 
the figure - the numerical integration results run parallel to the i -1 / 4 lines. The earliest time depletion zone 
behavior that precedes the Zeldovich regime, where pit) is expected to decay as t~ x / 2 , is more apparent in 
the system with the smaller diffusion coefficient. We have shown elsewhere that the transition from depletion 
zone behavior, t -1 / 2 , to Zeldovich behavior, t^ 1 / 4 , occurs at a time t s inversely proportional to the diffusion 
coefficient. The rapid decrease in the density due to the onset of finite size effects is clearly visible in both 
cases. For the larger diffusion coefficient (more rapid mixing) the onset of the classical regime is seen clearly, 
as indicated by the fact that the density decay curve in (a) runs parallel to the t^ 1 line. The time of the run 
in Fig. |l](b) is not sufficiently long to see this classical regime clearly. Note that, as expected, the finite size 
effects set in later for the system with the smaller diffusion coefficient. 

A point to note concerns the density moments of a strictly random mixture of ^4's and £Ts and those of 
a homogeneous mixture. In both of these the average of the product of the densities of different species is 
simply the product of the averages: (paPb) = (pa) (pb)- Consequently, an average of Eq. (|J) immediately 
leads to second order kinetics in both cases. It is interesting to note, however, that this equality of moments 
does not extend to other moments; in particular, the sum and difference variables for these two distributions 
behave quite differently. In a random distribution there are fluctuations in both the total (A plus B) local 
density and also in the difference variable, while in a homogeneous distribution these fluctuations have been 
smoothed out - indeed, this is what is meant by "homogeneous." For a random distribution the average 
of the square of the difference variable is linear in the average density, (7 2 (x)) = pmaxPi an d the average 
of the square of the sum variable also has such a linear contribution, cf. Eq. (|l3|) [Argyrakis et al., 1993; 
Lindenberg et al., 1996]. On the other hand, for a homogeneous distribution the difference variable vanishes 
so that (7 2 (a?)) = Oj while the average of the square of the sum variable is quadratic in the average density, 
cf. Eq. (|T4|). The kinetics in the reaction-diffusion picture depend only on the difference (p 2 ) — (j 2 ), which 
is the same for both cases. 

3.2 A detour 

The above description of the evolution of the reaction-diffusion system from an initial random distribution of 
reactants predicts the behavior of the global density in various kinetic regimes, Eq. (E^) followed by Eq. (p5|) 
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followed by Eq. (p7|), and in some cases predicts the approximate time around which a transition from one 
kinetic behavior to another takes place. However, we are not easily able to predict the actual crossover 
behavior itself. Thus, for example, we know when finite size effects set in, and we can reason that the global 
density will drop sharply around that time, but we have not predicted the functional form of the drop. 

Some insights can be obtained by pausing at this point to analyze the behavior of the different terms in 
the reaction-diffusion equation itself in some detail. This "generic" examination will help us understand not 
only the case of random initial conditions but also the behavior of the system with other initial conditions. 




The analysis is facilitated by considering the evolution of the global density from the following initial 
condition (more easily achievable than the random): A is initially placed in a single cell centered at xa and 
nowhere else, and B is initially placed in a different single cell centered at and nowhere else. In our 
(one-dimensional) numerical implementation of a system with N = 64 cells we place these two cells as far 
from one another as possible, i.e., there are 31 empty cells between the two initially occupied cells. The 
resulting evolution is shown in Fig. |[ The global density essentially does not decay at all until A's and £?'s 
on average encounter one another. In this numerical implementation the time for these encounters is roughly 
the same as the time at which finite size effects set in. There is an abrupt decrease in the density and the 
system then settles into classical behavior, as indicated by the straight t^ 1 line. This system is fairly trivial 
but offers an opportunity to make a number of observations about the various terms in the reaction-diffusion 
equation. 

The segregation function for this initial condition in which segregation is imposed from the outset is 
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easily calculated in the large L limit. In one dimension we find 



G{t) = -p= (l - e -(xA-x B ) 2 /8Dt^ (28) 



As long as t <C t m = (xa — xb) 2 /8D, the segregation function has the same time dependence, i -1 / 2 , as that 
of the random initial condition, viz. Eq. ( pl[ ) let us for the moment concentrate on this time regime. As in 
the random case, the decay of the driver G(t) is sufficiently slow to be a dominant contribution in Eq. (|ll|). 
As in the random case, this dominant contribution must be balanced by the other dominant contribution, 
(p 2 ). However there is a major difference between this case and the random one: in the case of a random 
initial condition (p 2 ) is related in a simple fashion to the global density p(t), cf. Eqs. ( |l3| ) and (jlj). Therefore 
the balance of the two terms on the right of Eq. (|ll|) immediately tells us the kinetic behavior of the global 
density, as described earlier. That led to the depletion zone behavior at early times (i -1 / 2 ) followed by the 
Zeldovich behavior at later times (t -1 / 4 ). The balance between p(i) and the difference between (p 2 ) and 
(j 2 ) simply confirms this kinetic behavior. 

In the present example there are two (related) differences from the random initial condition case. First, 
there is no simple relation between the average of the square of the sum variable and the global density. 
Indeed, for times t <C t m the reaction plays essentially no role, so that one can quite accurately evaluate 
(p 2 ) just from the diffusion portion of the reaction-diffusion equation to obtain 



(p 2 (x, t)) » -^L ( 1 + e -(^-« ) 2 / 8D *) . (29) 
\l Dt ^ 



On the other hand, the global density is constant in this time regime, so that there is no simple moment 
relation between (p 2 ) and the global density. The leading contributions of Eqs. (|2^) and ( |29| ) cancel as they 
do in the random case, but the difference between them is now exponentially small. Here this in itself says 
nothing about the decay of the global density, since this difference is not related in any simple way to the 
global density. Information about the global density is obtained only when we explicitly balance p(t) against 
this very small difference. 

In our numerical example the time t m when reactants encounter one another on average is of the same 
order as the time tf. At around this time the segregation function G(t) drops abruptly. So does the global 
density, which then settles to the classical decay law (iff) -1 . When this happens, the segregation function 
becomes irrelevant, the moment relation Eq. (|l4|) becomes applicable as diffusion homogenizes the system, 
and the principal balance in Eq. ([ll]) is between p(t) and the first term on the right, that is, Eq. (|ll|) reduces 
to the classical rate law. 

This description together with that of the random initial condition case permits us to reach some generic 
conclusions about the behavior of the different contributions to Eq. ([TT|). Figure || is a schematic (not to 
any particular scale) of the evolution of these different contributions in low dimensions. First consider the 
segregation function G(t) = (j 2 ^ (solid curve). It typically begins by decaying according to a power law, 
more slowly than t -2 , with perhaps a transition from one slow power decay to another. At some time, 
at the latest the time tf when finite system size effects set in, the segregation function drops abruptly 
(exponentially) as the system becomes homogenized. Next consider the average of the square of the density 
(p 2 ) (dashed curve). It is larger than the segregation function but decays to leading order in the same fashion 
before time tf. These two contributions are the leading terms in Eq. ( |Tl| ) and therefore balance one another. 
The difference between them decays more rapidly than either term, so in this regime p(t) (which balances 
this difference and is sketched as a dot-dashed curve) decays rapidly as well and is negligible in the leading 
balance. The square of the global density, p 2 (t), is also sketched (dotted curve). It is generally smaller than 
the other two contributions, but approaches (p 2 ) (for different reasons in different cases). As it does so, the 
global density itself is determined by the leading balance, that is, p(t) ~ G 1 / 2 (i), and this relation must of 
course be consistent with the balance between p and the difference between the leading terms. During the 
crossover time when G(t) drops abruptly the system homogenizes, and p catches up with p 2 , which by then 
is essentially equal to (p 2 )- G(t) then becomes irrelevant and the kinetics is classical. 

This generic sequence describes every initial condition that we have dealt with here and elsewhere [Ar- 
gyrakis et al., 1993; Lindenberg et al., 1994, 1996; Sancho et al., 1996]. The differences arise in the specific 



11 



10° 




100 1000 10000 

t 



Figure 3: Schematic of the typical behavior of the various contributions to Eq. ([□]): G(t) = (~/ 2 ), (p 2 ), p 2 , 
-p/K. 



slopes associated with a variety of anomalous kinetics (that is, decays of the global density that are slower 
than the classical t~ 1 ) before finite size effects set in. These anomalous kinetics arise from (indeed, require) 
initial density fluctuations and show a range of behavior depending on the distribution of these fluctuations. 
The random initial condition is but one example, and in one dimension is associated with decays of the global 
density of the form t~ x / 2 followed by < -1 / 4 . In the examples to follow we deal with the "more anomalous 
case" of a fractal initial distribution leading to a slower decay than that of the random case, and with a "less 
anomalous" case of initially correlated pairs leading to a decay closer to the classical. 

3.3 Fractal initial distribution 

Consider next an initial condition in which the long wavelength components of the difference variable are 
emphasized relative to the random initial distribution, leading us to expect a greater deviation from classical 
behavior [Provata et ai, 1996]. 

Our initial distribution is generated as follows: the number N of cells in a row is chosen to be a power 
of 3, N = 3". We divide N by 3, thereby generating JV/3 domains (each containing many cells). One of 
the three domains, chosen at random from among the three, remains empty and is discarded from further 
consideration. The other two domains are each subdivided into three subdomains, and one of each group 
of three subdomains again remains empty and is set aside, while the other two are again subdivided. The 
process continues in this way until the selected domains are as small as possible, that is, of size Ax. An 
initial density po of one of the reactants is placed in each of the selected elementary cells. The process is 
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repeated for the second reactant, and these two distributions are superimposed. In the limit N — > oo this 
leads to a fractal geometry of occupied cells of dimension Df = In 2/ In 3 = 0.63 [Vicsek, 1992]. The variable 
7 is also characterized by the same initial fractal dimension, and the structure function for small k behaves 
as a power law with an exponent given by the fractal dimension, 

S(k) ~ k-°f (30) 

Note that the exponent is negative, i.e., the long wavelength components are enhanced relative to those of a 
random intitial distribution. 

Consider now the sequence of behaviors that might be expected with this initial condition. Since there 
may be some premixed cells, there will be a rapid decay of density in these cells on a time scale K . Once 
this is over, we expect a lull in the activity since aside from the premixed cells, occupied cells are in general 
quite separate. Indeed, if the fractal dimension is too high (e.g. if instead of subdividing into 3 parts at each 
generation we subdivided into m portions with m fairly large) then there are many premixed cells and much 
of the reactant disappears during the very early time regime. The quiescent period then dominates the time 
evolution because there is little reactant left and it takes a very long time for reactants to find one another 
(somewhat like our point source example). 

During the quiescent period the distribution is rather segregated but the sum variable is becoming 
spatially homogenized through diffusion; hence one expects the moment relation (|l4| ) to apply. The difference 
variable, on the other hand, is determined by the initial fractal pattern and leads to the segregation function 
determined by Eq. ( |30|) . The rate is then of the form 

^ ~ -K[p\t) - Q(Df)- (d ~ D/)/2 ]- (31) 

For d = 1 (and, indeed, for d < 4 + Df) the dominant terms are those on the right hand side, and the global 
density thus decays as 

p((t) ~t- {1 - D f^ 4 . (32) 

Note that this decay is even slower than the Zeldovich behavior. This decay continues until finite size effects 
and the subsequent behavior described earlier take over. 

In Fig. H w e show the results of the numerical integration of the reaction-diffusion equations and the 
prediction (|3^). It should also be noted that our initial patterns are not "ideal" fractals because we stop 
the subdivision at a finite cell size. Our systems have N = 81, so we only go through four generations of 
subdivisions. Nevertheless, the predicted results are clearly apparent in the figures. In Fig. [|(a) we see the 
slow decay of the form t~ 092 expected for an initial distribution of fractal dimension Df — 0.63, followed by 
the abrupt decrease in the global density associated with finite system size effects and the eventual approach 
to classical behavior. Figure [|(b) is for a smaller diffusion coefficient - the slow decay continues for a longer 
time, finite system size effects set in later, and the time scale is not sufficiently long to see the eventual 
appearance of classical behavior. 



3.4 Initially correlated pairs 

A generic example of a constrained initial distribution of reactants in which the long wavelength components 
of the initial density difference are diminished through built-in correlations is that in which randomly oriented 
A — B pairs of particles separated by a fixed distance c are initially distributed randomly in the system. 
Subsequent to this initial constrained placement, each particle moves freely. Lindenberg et al. [Lindenberg 
et al., 1994, 1996; Sancho et al., 1996] have studied in detail the kinetic progression for the A + B — > 
reaction with this initial distribution. In our discretized reaction-diffusion picture we translate this into the 
initial random placement of pairs of cells with centers separated by a distance A (c = A Ax), one cell of each 
pair containing only species A and the other only species B (whether or not we allow initial occupancy of a 
cell by different species arising from different pairs is immaterial - it only changes the very early transient, 
as we have noted above). 

The segregation function for this initial condition in an infinite volume is given by [Lindenberg et al., 
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Figure 4: Kinetics for the fractal-like initial pattern described in the text. In both figures N = 81 and 
K = 10. (a) D — 2.0, (b) D = 0.5. The curves are the numerical solutions of the reaction-diffusion equations. 
Straight lines in (a): t^ 1 - D f^ i = t' - 092 and t' 1 . Straight line in (b): t~ a - 092 . 



1994, 1996] 

G(t) = Q{Dt)- d ' 2 (l - e - x2 / SDt ^j , (33) 

where, as before, Q is a constant proportional to the initial density. At times t < t c = A 2 /8D the segregation 
function behaves exactly as in the case of the random initial distribution. Beyond time t c the decay of G(t) 
changes from t~ x l 2 to t~ 3 / 2 (the associated behavior of the structure function for small k is S(k) ~ k 2 ) 
[Sancho et at, 1996]. This small k behavior reflects the constraints on the long wavelength components of 
the initial density difference fluctuations imposed by the pair correlations. 

The decay of the global density depends, as always, on the balance of terms in Eq. ([ll]) and on the 
relation between (p 2 ) and the global density. 

In order to describe the ensuing kinetic regimes, we must distinguish two cases depending on the relation 
between the time t s and the time t c . Recall from our discussion of the random initial condition that t s is the 
crossover time at which the relation between the second and first moment of an initially random distribution 
changes from the form Eq. ( |l3| ) to Eq. (|lj), that is, for the segregation process to begin. If t c < t s then 
segregation never sets in because initially correlated pairs of cells react (in a "burst") before segregation 
can begin. We have discussed this case in detail in [Lindenberg et ai, 1994, 1996] and will not present this 
particular case here. If t c > t s then up to time t c the system evolves as would a random distribution, with an 
early p ~ t -1 / 2 depletion zone evolution followed by a p ~ i" 1 / 4 Zeldovich evolution during which segregated 
regions grow in size. These two regimes are described by Eqs. (|2^) an d (pi[ ) respectively. However, at around 
time t c the kinetic behavior changes because the second term in Eq. (|33|) becomes important: as t increases 
well beyond i c , instead of Eq. ( p4| ) we now have for d < 2 a rate law of the form 

^ ~ ~K[p 2 {t)~Q{Dt)-( d+2 V 2 ]. (34) 
at 
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The global density in one dimension now decays as 



Piit) 



t 



-3/4 



(35) 



We have explained the origin of this behavior (note that this decay, although anomalous, is faster than that 
of an A + A reaction): the distance A effectively sets a limit on the sizes of the segregated aggregates. The 
initial correlations in effect homogenize the system by constraining the fluctuations of wavelength greater 
than A. The relatively rapid decay Eq. (|3^) occurs because after time t c the single species aggregates actually 
begin to shrink, ft is this shrinkage that is reflected in the more rapid decay. 

Eventually system size effects enter the picture, Eq. (|3^) can no longer be used, and G(t) decays, as 
always, exponentially. The kinetic evolution beyond this point is classical. 
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Figure 5: Kinetics for initially correlated pairs. N = 4096, K = 10, D = 2, and A = 32. The curve is the 
numerical solution of the reaction-diffusion equations. Straight lines: i -1 / 4 and i -3 / 4 . 



Figure || illustrates the salient features of this behavior. We first see a decay of the global density close 
to the Zeldovich i -1 / 4 , followed by a crossover to a i~ 3 / 4 aggregate shrinkage decay beginning at around the 
time t c ~ A 2 /8D = 64. Note that for the parameters used in this figure tf is off the scale and finite size 
effects are thus not seen. 
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4 Conclusions 



We have presented a study of the effects of initial distributions on the behavior of the irreversible reaction 
A+B — ► under stoichiometric conditions, mainly in one dimension (d = 1). Our work is based on a reaction- 
diffusion model, and we include detailed discussions of the behavior of each term in the reaction-diffusion 
equations. We have focused on the effects of the initial pattern of fluctuations of the density difference and 
have stressed that this pattern determines the evolution of the system for all time. In particular, the long- 
wavelength components of the initial difference distribution determine the asymptotic decay of the reactant 
densities. Thus, this system constitutes a temporal "mirror" that for all time reflects the instantaneous 
initial spatial pattern. 

An initially random distribution of reactants leads to well-known and amply studied segregation of species 
and the associated slowing of the reaction. In an infinite system the asymptotic decay of the reactants in this 
case is described by the "Zeldovich law" p <~ i~ d / 4 = i -1 / 4 , which is distinctly different from the classical 
law-of-mass-action behavior p <~ t -1 . Initial distributions that limit the long wavelength components of 
the difference fluctuations lead to constraints on the macroscopic segregation effects and to rate laws closer 
to the classical. We reviewed a particular example of such an initial distribution, namely, one containing 
initially correlated A — B pairs. On the other hand, initial distributions that enhance the long wavelength 
components lead to even greater deviations from the law of mass action. For instance, an initial fractal 
difference pattern of dimension Df leads to the very slow de cay p ~ i-( d - D /)/ 4 . 

We have emphasized the fact, also well known, that no single exponent characterizes the density decay 
for all time, even in an infinite system. Instead, we have indicated that it is often possible to characterize 
the decay by a single exponent for a considerable length of time (compared to the other time scales in the 
problem as determined by the reaction rate coefficient and the diffusion coefficient), followed by a crossover 
of relatively short duration that leads to another exponent for a considerable time span, etc. The sort of 
analysis presented here is only appropriate when these time scale separations are possible. 

Our numerical work is based on a direct integration of the reaction-diffusion equations and is, in all cases, 
subject to the long wavelength constraints imposed by a finite system size. Thus, although the anomalous 
decays indicated above might be expected to persist forever for irreversible diffusion-limited reactions in 
infinite systems, in the present work we do not emphasize the idealized asymptotic behavior. Instead, we 
include a discussion of the behavior beyond the time tf at which finite system size effects become apparent. 
The reaction-diffusion model predicts classical behavior at long times. We should note that although in the 
paper we discuss averages over initial distributions, the numerical results presented in our figures are all 
for single realizations of these distributions. Further averaging over initial distributions did not modify the 
outcomes. 

It would be interesting to observe in detail not only the time evolution of the reactant densities but also 
the spatial patterns associated with this evolution for the different initial distributions that we have studied. 
We plan to do so in the future. 
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